### Alizade, Dancygier, Ditlmann 
### "National Penalties Reversed"
### Replication Code 
### Table 3
### For questions, contact jalizade@princeton.edu

# empty environment
rm(list = ls())

#setwd("")
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")

# load necessary packages
library(readstata13)
library(boot)

# load data set
dat <- read.dta13("data_experimental.dta")


### table ###

# treatment
treat <- rep(c("Nationality", "Citizenship", "Motivation"), each = 2)


# treatment condition
treat_cond <- c("Canadian", "Turkish", "Single", "Dual", "Feel German", "Legal Benefits")


# N
n <- c(sum(dat$e1_treat_turkish==0, na.rm = T), sum(dat$e1_treat_turkish==1, na.rm = T),
       sum(dat$e1_treat_dual==0, na.rm = T), sum(dat$e1_treat_dual==1, na.rm = T),
       sum(dat$e1_treat_benefits==0, na.rm = T), sum(dat$e1_treat_benefits==1, na.rm = T))


# response rate
resp_rate <- c(mean(dat$e1_response[dat$e1_treat_turkish==0], na.rm = T), mean(dat$e1_response[dat$e1_treat_turkish==1], na.rm = T),
               mean(dat$e1_response[dat$e1_treat_dual==0], na.rm = T), mean(dat$e1_response[dat$e1_treat_dual==1], na.rm = T),
               mean(dat$e1_response[dat$e1_treat_benefits==0], na.rm = T), mean(dat$e1_response[dat$e1_treat_benefits==1], na.rm = T))
resp_rate <- round(resp_rate*100, 2)


# ratio
ratio <- c(mean(dat$e1_response[dat$e1_treat_turkish==1], na.rm = T) / mean(dat$e1_response[dat$e1_treat_turkish==0], na.rm = T), NA,
           mean(dat$e1_response[dat$e1_treat_dual==1], na.rm = T) / mean(dat$e1_response[dat$e1_treat_dual==0], na.rm = T), NA,
           mean(dat$e1_response[dat$e1_treat_benefits==1], na.rm = T) / mean(dat$e1_response[dat$e1_treat_benefits==0], na.rm = T), NA)
ratio <- round(ratio, 2)


# bootstrap CI for ratio
set.seed(9447)

ratio_turk <- function(dat, i) {
  dat2 <- dat[i,]
   mean(dat2$e1_response[dat2$e1_treat_turkish==1], na.rm = T)/mean(dat2$e1_response[dat2$e1_treat_turkish==0], na.rm = T)
}
boot_ratio_turk <- boot(data=dat, statistic=ratio_turk, R=10000)
boot.ci(boot_ratio_turk, conf = 0.95, type="perc")
ratio_ci_turk <- round(boot.ci(boot_ratio_turk, conf = 0.95, type="perc")$percent[4:5], 2)

ratio_dual <- function(dat, i) {
  dat2 <- dat[i,]
   mean(dat2$e1_response[dat2$e1_treat_dual==1], na.rm = T)/mean(dat2$e1_response[dat2$e1_treat_dual==0], na.rm = T)
}
boot_ratio_dual <- boot(data=dat, statistic=ratio_dual, R=10000)
ratio_ci_dual <- round(boot.ci(boot_ratio_dual, conf = 0.95, type="perc")$percent[4:5], 2)

ratio_benefits <- function(dat, i) {
  dat2 <- dat[i,]
   mean(dat2$e1_response[dat2$e1_treat_benefits==1], na.rm = T)/mean(dat2$e1_response[dat2$e1_treat_benefits==0], na.rm = T)
}
boot_ratio_benefits <- boot(data=dat, statistic=ratio_benefits, R=10000)
ratio_ci_benefits <- round(boot.ci(boot_ratio_benefits, conf = 0.95, type="perc")$percent[4:5], 2)

# combine CIs into one data frame
ratio_ci <- rbind.data.frame(ratio_ci_turk, c(NA, NA), ratio_ci_dual, c(NA, NA), ratio_ci_benefits, c(NA, NA))
names(ratio_ci) <- c("ratio_CI_lower", "ratio_CI_upper")


# difference in means
diff <- c(mean(dat$e1_response[dat$e1_treat_turkish==1], na.rm = T) - mean(dat$e1_response[dat$e1_treat_turkish==0], na.rm = T), NA,
          mean(dat$e1_response[dat$e1_treat_dual==1], na.rm = T) - mean(dat$e1_response[dat$e1_treat_dual==0], na.rm = T), NA,
          mean(dat$e1_response[dat$e1_treat_benefits==1], na.rm = T) - mean(dat$e1_response[dat$e1_treat_benefits==0], na.rm = T), NA)
diff <- round(diff*100, 2)


## standard error of difference
# perform t-tests first
diff_t <- c(t.test(dat$e1_response ~ dat$e1_treat_turkish)$statistic, NA,
            t.test(dat$e1_response ~ dat$e1_treat_dual)$statistic, NA,
            t.test(dat$e1_response ~ dat$e1_treat_benefits)$statistic, NA)

# to obtain the standard error, divide the mean difference by the t-statistic
diff_se <- round(abs(diff / diff_t), 3)


## p-value of difference
# perform t-tests first and store the p-values
diff_p <- c(t.test(dat$e1_response ~ dat$e1_treat_turkish)$p.value, NA,
            t.test(dat$e1_response ~ dat$e1_treat_dual)$p.value, NA,
            t.test(dat$e1_response ~ dat$e1_treat_benefits)$p.value, NA)
diff_p <- round(diff_p, 3)


# generate table
table <- data.frame(treat, treat_cond, n, resp_rate, ratio, ratio_ci, diff, diff_se, diff_p)
table


